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The sign problem is the fundamental limitation to quantum Monte Carlo simulations of the 
statistical mechanics of interacting fermions. Determinant quantum Monte Carlo (DQMC) is one 
of the leading methods to study lattice models such as the Hubbard Hamiltonian, which describe 
strongly correlated phenomena including magnetism, metal-insulator transitions, and (possibly) 
exotic superconductivity. Here, we provide a comprehensive dataset on the geometry dependence 
of the DQMC sign problem for different densities, interaction strengths, temperatures, and spatial 
lattice sizes. We supplement these data with several observations concerning general trends in 
the data, including the dependence on spatial volume and how this can be probed by examining 
decoupled clusters, the scaling of the sign in the vicinity of a particle-hole symmetric point, and the 
correlation between the total sign and the signs for the individual spin species. 


PACS numbers: 71.10.Fd, 02.70.Uu 

I. INTRODUCTION 

Monte Carlo simulations of classical systems have 
generated very precise information about phases and 
phase transitions in statistical mechanics. One dramatic 
example of the power of the methodology is that of the 
Ising model, where the transition temperature on a cubic 
lattice is now known to six decimal places^ and critical 
exponents have been evaluated to four decimal placesi^ A 
roughly similar situation holds for unfrustrated quantum 
spin and boson systems. For example, in the Bose- 
Hubbard Hamiltonian^^ the critical interaction strengths 
at a density of one boson per site for the superfluid 
to Mott insulator transition in the ground state, are 
available^ to an accuracy of better than one part in 
10^. Spatial lattice sizes are somewhat smaller than for 
classical problems since writing the partition function as 
a path integral introduces an additional ‘imaginary time’ 
dimension, but are nevertheless quite large, e.g., up to 
10"^ sites for the Bose-Hubbard example cited above. 

Fermions (in more than one dimension) are more 
challenging for two reasons. First, the fermionic action 
is non-local: the Boltzmann weight typically takes the 
form of a determinant. Thus, updating all the degrees 
of freedom has a computation time which scales as 
a nonlinear power of the system size N. The cost 
of a method like determinant quantum Monte Carlo 
(DQMC)^ scales as N^, as opposed to the naive linear 
in N scaling (ignoring such complications as critical 
slowing down) in many classical and quantum spin/boson 
applications. Second, and far worse, there is no guarantee 
the sign of the determinant, which is used as the 
probability, is positive. Although one can formally use 
the absolute value as a weight, and include the sign in 
the measurement, in practice one ends up evaluating 
the ratio of two numbers, which becomes dominated by 
statistical error at low temperatures as they both become 
very small. This situation is known as the “fermion sign 
problem” whose solution is conjectured to be NP- 
hardi^ At present, therefore, there is no known method of 


accessing the low temperature properties of Hamiltonians 
like the Fermi-Hubbard model with QMC. 

Fortunately, there are some situations where the sign 
problem is not manifest. For example, the Boltzmann 
weight often takes the form of the product of two 
determinants, one for each of the two electron spin 
species, and it can happen that the signs of the individual 
determinants perfectly match, so that the Boltzmann 
weight remains positive. This occurs in the complete 
parameter range of the Hubbard model with attractive 
interactions, enabling a study of superconductivity and 
charge density wave physics. It also occurs in the 
Hubbard model with repulsive interactions on bipartite 
lattices in the limit of average one particle per lattice 
site (half-filling), so that the Mott transition and long- 
range antiferromagnetism can be explored. Other 
instances of situations where the sign problem is absent 
are mentioned in Sec. HC. However, many of the 
most interesting questions concerning strong correlation 
physics remain inaccessible, most notably the question 
of whether the two-dimensional (2D) square lattice 
repulsive Hubbard Hamiltonian has a low temperature d- 
wave superconducting transition, so that it would provide 
a good description of cuprate superconductivityi^ 

We have two main goals in this manuscript. The 
first is to present a set of data for the sign problem in 
DQMC for different geometries. These include bipartite 
[one-dimensional chain, ladder, 2D square, and three- 
dimensional (3D) cubic, Lieb, honeycomb and the 1/5- 
depleted square] lattices, and non-bipartite Kagome 
and triangular lattices. We consider DQMC because 
it is a powerful and widely utilized approach to the 
correlated electron problem. Our second goal is to 
discern trends in the DQMC sign problem. We will 
consider, for example, several new issues: how the sign 
depends on the spatial lattice size, the scaling of the 
sign in the vicinity of particle-hole symmetric (PHS) 
points, and the ‘entanglement’ of the sign as probed 
by the consideration of coupled cluster geometries. 
Although statements concerning the first of these points. 
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the spatial size dependence, have been made in the 
literature, numerical data are rather scanty, owing to 
the computational limitations existing in the initial 
investigations. Specifically, the data in Refs. 0 and M 
were restricted to 4x4, 6x6, and 8x8 lattices. As a 
consequence, the scaling regime was not reached for many 
of the parameter sets. For example, the average sign 
sometimes increases as the system size grows, rather 
than decreasing. This situation is rectified here. The 
lattices we study also have various unique features, such 
as nesting of the Fermi surface, van-Hove singularities 
in the density of states N{oj) = J2k ~ and flat 
bands, whose possible effects on the sign problem we will 
examine. 


It is worth noting what we will not cover: There 
are a number of methods which are closely related to 
DQMC in that they involve a Hubbard-Stratonovich 
decoupling of the interaction, and a Boltzmann weight 
built from fermion determinants. These include impurity 
algorithms as well as the dynamical mean field 
theory (DMFT)i^“— and its cluster extensions, the 
dynamical cluster approximation (DCA)fi^ and the 
cellular DMFTii^ A strength of some of these methods 
is that the sign problem is greatly mitigated relative 
to DQMC, at least if the cluster size is not too large. 
This is true even if the bath degrees of freedom are 
discretizedjii Also closely related to DQMC are zero- 
temperature algorithms which use e~^^ as a projection 
operator on a trial wave functionJ^d^ Constraints can be 
introduced in these ground state methods to eliminate 
the sign problem, at the expense of systematic errors 
in the solutionii^ Despite their relations to DQMC, we 
will not discuss these approaches here. Similarly, within 
DQMC itself there are different choices of the manner in 
which the Hubbard-Stratonovich field is decoupledi^^— 
Here, we base our calculations on only the “density 
decoupling”, described in Sec. HB. Though we do 
not explicitly consider the above related approaches, we 
expect that some of our results and general analysis may 
have applications to them as well. 


The remainder of this paper is organized as follows: 
In Sec. H we review the Hubbard Hamiltonian and the 
basic formulation of DQMC, followed by some general, 
well-known properties of the resulting sign problem. In 
Sec. HI we record values for the average sign for different 
lattice geometries. Section IV examines general patterns 
in this data. Finally, Sec. V contains concluding remarks. 


II. GENERAL CONSIDERATIONS 
CONCERNING THE SIGN PROBLEM 

A. The Hubbard Hamiltonian 

Our focus is on the single-band Hubbard Hamiltonian, 

H = -^Uo{ cL 

(ij)a io- 

+ ( 1 ) 

Here cJ^(Cj^) is the creation (destruction) operator for 
a fermion with spin a on site i = 1,2, - ■■ N and 
is the number operator, tij is the hopping 
amplitude between nearest-neighbor sites i and j, U is 
the interaction strength, and n is the chemical potential. 
For geometries where there is only one type of hopping, 
we set tij = t = 1 as the unit of energy. We will 
denote the first line of the Hamiltonian H in Eq. (HD 
by K, and the second line by V. This latter term is 
written in PHS form. (See Sec. HC.) The different models 
considered in this paper are distinguished solely by the 
geometry encoded in the near-neighbor designation (ij) 
in the kinetic energy term. Even with a common choice 
t = 1, different geometries have distinct bandwidths W, 
the spread of eigenvalues of the U = 0 (single particle) 
Hamiltonian. Although it is sometimes the case^^ that 
using W as the scale of kinetic energy, rather than t, 
produces better comparisons across different models, we 
did not find that to be useful here. We retain the 
standard convention of normalizing to t. 

There is, of course, much interest in generalizations 
of the Hubbard Hamiltonian, e.g., to multiple bands, 
longer-range density-density interactions, and Hund’s 
rule type interactions. Multiple bands, can in fact be 
written in the form of Eq. (HD, with the understanding 
that the label i incorporates both spatial and band 
indices. Thus, from the viewpoint of a DQMC 
simulation, setting up the 2D ‘periodic Anderson model’ 
(PAM) which has a square lattice of spatial sites 
and two orbitals per site, is formally identical to a 
two layer geometry, in which there is a single orbital 
on each site. Thus, with the freedom to choose 
Ui (tij) to be site/orbital (bond) dependent, Eq. ([T]) 
incorporates Hamiltonians like the PAM. Concerning 
intersite (interorbital) and Hund’s rule interactions, the 
sign problem is typically much worse than for an on¬ 
site U between fermions of different spin species only. 
For example, for a model of the Cu02 planes of 
cuprate superconductorsit was found that the sign 
problem restricted simulations to interactions Upd 1 at 
temperatures where local spin correlations were seen to 
begin to develop. Like DQMC, Hund’s rule interactions 
also present grave sign problem difficulties in DMFTj^ 
We do not explicitly consider them here. 
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B. Determinant Quantum Monte Carlo 


The fundamental idea of DQMC— is to take 
advantage of the fact that it is possible to compute 
analytically the trace of a product of the exponentials 
of quadratic forms of fermion creation and destruction 
operators. If we denote the vector of creation operators 
(c|, C 2 , Cg • • • c]y ) by c, and Aj are (symmetric) N x N 
matrices of real numbers, then 

Tr ( gcMac . . . gcMi,c ^ 

= det{I + BiB2B3---BL) ( 2 ) 

Here Bi = It is important to emphasize that 

the trace in the left hand side of Eq. m is over a 
2 ^ dimensional Hilbert space of the fermionic operators 
while the determinant on the right hand side is taken 
over a real matrix of dimension N. 

The interaction term in the Hubbard Hamiltonian 
Eq. (m is not quadratic in the fermionic operators, 
but can be made so by first discretizing the inverse 
temperature /3 = Lt and then employing the Trotter 
approximation^l^— 

Tr = Tr ( ) (3) 

RiTr(e e e e •••e e ) 


Here the exponential of the full Hamiltonian H = K + V 
is approximated by the product of the exponentials of K 
and E, a well-controlled procedure which can be made 
arbitrarily accurate by taking r —^ 0 . 

The purpose of this procedure is the^ isolation of 
the exponential of the interaction term V, which can 
then be rewritten using a Hubbard-Stratonovich(HS) 
transformation: 

g-r£/(nit-5)(nij.-|) _ lg-f7T/4 ^ gA A" (ni^-^ 

^ A'=±l 


where cosh A = Because one needs to transform 

the interaction term on every spatial site i = 1,2, ■■■ JV 
and also for each of the I = 1,2, - ■■ L exponentials of tV 
in Eq. there are a total oi NL HS variables 

In Eq. (IH), we have employed the discrete HS 
transformation introduced by Hirsch;^ but one could 
also use a continuous variable X and a Gaussian integral. 


g-rC/(nit-5)(ni4.-|) 





(5) 


with 7 = ^Ut 12. There are some differences in the 
efficiency of the exploration of phase space between the 
discrete and continuous cases.— 

Once the HS transformation is introduced, all the 
exponentials in the trace of Eq. (|31) are quadratic in the 
fermion operators, so, the identity in Eq. @ can be used 


to perform the trace over the Hilbert space analytically. 
The sum over the HS configurations X{i,l) is performed 
stochastically using Monte Carlo techniques. The 
corresponding Boltzmann weight takes the form of the 
product of two determinants (one for each spin specie) of 
matrices M^{X) of dimension N. As this determinant 
product may be negative for some HS configurations, 
the sampling is done using the absolute values of the 
determinant product, and measured expectation values 
are adjusted accordingly. 

The average sign (S) is then defined to be the ratio 
of the integral of the product of up and down spin 
determinants, to the integral of the absolute value of the 
product. An analogous definition holds for the average 
sign {Sa) of the individual determinants: 


^ Ea:' detMi-(A’)detM;(A’) 

^ Ex |detA/t(A’)detM^(A)| 
/ q \ _ Ex detAfo.(A’) 

^ Ea |detM<,(A)r 


( 6 ) 


In the case we consider here, with no external magnetic 
field, by symmetry (S-^) = (Si). 

As a practical matter, these quantities are obtained 
by generating configurations with the (non-negative) 
weight |detM^(A’)detM^(A’)| and measuring the ratios 
detM^(A’) detM 4 ,(A’) / | detM^(A’) detM 4 ,(A’) | and 

detMcr(A’) / |detMCT(A’)| for each configuration X. In 
Eqs. o and ®, we have coupled the HS variable to 
the 2 ; component of fermionic spin, rii-i- — rii^. It is 
also possible to write transformations which involve the 
xy components of a spin, or even local pairing 

operators, These in general worsen the sign 

problem— In the attractive Hubbard Hamiltonian, 
the HS variable couples to the charge ni-\ + riii 
site i. This makes the matrices Ai identical for up 
and down fermions, so that their determinants are also 
identical, and thus, there is no sign problem. If a charge 
decoupling is used for the repulsive model, the HS 
transformation would involve complex numbers, and the 
determinants would be complex as well, leading to an 
even more challenging ‘phase problem’. 

We note that there are many details omitted in this 
brief description, including methods to stabilize the 
product of the Bi matrices so that round-off errors 
do not accumulate, the precise Monte Carlo update 
procedure (how many variables are altered in each step), 
more rapid procedures for obtaining the ratio of new 
to old determinants after a HS variable is updated, 
how to evaluate non-equal time observables, analytic 
continuation to obtain dynamic behavior, and so forth. 
The reader is referred to Refs. for more 

complete discussions. 
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C. Particle-Hole Symmetry 

On a bipartite lattice, and at ^ = 0, the Hamiltonian 
Eq. ([T]) is PHS. That is, the Hamiltonian is invariant 
under the transformation —>■ (—l)®Cjg.. Here (—1)* = 
+1(—1) on the A(B) sublattice. This symmetry is present 
even if t and U vary spatially across the lattice. As one 
physical consequence, the density pi = (+ riii ) = 1 
(half-filling) for all values of t, U and T. Correlation 
functions at density p and 2 — p have the same values, 
or are trivially related. Particle-hole transformations 
involving only one spin specie can also be used to relate 
the attractive and repulsive Hubbard Hamiltonians in 
this limit. 

PHS has profound implications for DQMC. When 
it is present, the determinants of the up and down 
matrices have the same sign. Thus, although they 
individually can go negative, their product is always 
positive. As a consequence, low temperature physics 
can be accessed at half-filling. This fact enabled DQMC 
to establish rigorously^^iM that the single band, square 
lattice Hubbard model has long-range Neel order at 
T = 0, as opposed to a disordered (resonating valence 
bond) ground state. 

The sign problem can be absent in some other types of 
Hamiltonians with a similar symmetry requirement, for 
example in a model with an interaction which takes the 
form of the square of near-neighbor hopping^ for a low 
energy theory of the onset of antiferromagnetism,— and 
in spin polarized Fermi systems.— Indeed, the number 
of special situations where the sign problem is absent 
is rapidly growing, including, for example, in quasi-lD 
condensed matter models of ferromagnetism;^ and via 
the ‘fermion bag’ approach, in lattice gauge theor y 
The sign problem is also absent in Hubbard models 
with a larger number of spin components,— and, very 
interestingly, in a class of spinless fermion models^, a 
unique situation where positivity is not dependent on 
having an even number of fermionic species. 

III. SIGN PROBLEM DATASETS 

Our goal in this section is to present a unified and easily 
comparable collection of data for the sign problem for 
different lattice geometries, including hypercubic lattices 
in dimensions d = 1,2,3; other bipartite structures like 
the honeycomb, Lieb and 1/5-depleted square lattices; 
and finally two non-bipartite lattices: triangular and 
Kagome. For each case we will exhibit the average sign 
for a range of temperatures T, interaction strengths U, 
lattice sizes, and densities. 

We focus on the product of the signs of the up and 
down determinants, since that is what is relevant for 
extracting physics from the DQMC simulation. However, 
in Sec. IVB we will present a brief analysis of the 
individual spin components. Among other things, we 
observe that, except at PHS points, the signs of the up 


and down determinants tend to be rather uncorrelated, 
so that the signs of the individual components lend no 
further information, and can be qualitatively inferred 
from the square root of the total sign. 

Some of these geometries have unique features in their 
non-interacting densities of states (DOS). The square 
lattice possesses a van-Hove DOS singularity at p = 1. 
In contrast, the honeycomb lattice has a DOS which 
vanishes linearly there. The Lieb lattice has a flat energy 
band between two dispersing ones, while the flat band in 
the Kagome lattice can be chosen to be either the lowest 
or highest set of energy levels, depending on the sign of t. 
One of our goals is to examine how such features might 
affect the fermion sign, an issue to which we will return 
in the conclusion. 


A. Hypercubic Lattices 

In this subsection we present data for hypercubic 
lattices; linear chains, ladders, the square lattice, and 
the cubic lattice. In all cases, we use periodic boundary 
conditions except for the rungs of the ladder geometry. 

Although the sign does not cause a problem for 
Hubbard world line method s in one dimension, 
DQMC does have a sign problem in this case^. In Fig. [U 
the average fermion sign, {S) (also denoted (sign)), for 
the chain geometry is shown for fixed {7 = 4 and /3 = 8. 
Figure [T](a) shows the doping dependence of the sign, 
which has a non-trivial structure. Most notably, it shows 
a peak around p = 0.5. One may wonder if such a local 
maximum arises due to remnants of the ‘shell’ effect. 
That is, at {7 = 0 the k space grid is sufficiently coarse 
such that p{p) shows distinct plateaus where (S) tends 
to be closer to one. This phenomenon is well known, 
for example, on square lattices that are not too large 
(e.g., 4 X 4). However, as will be seen in Fig. O this 
occurs only on small lattices, and is unlikely the origin 
of the maximum at quarter filling here where the k space 
grid is much more refined. Another interesting feature 
is that (S') remains small at low density. This is, again, 
rather different from what happens on a square lattice 
where (S) —>■ 1 as p —>■ 0 (Fig. or even ladders (Fig. [2]). 
The low value of (S) as p —0 does not appear to be 
connected to the divergence of the density of states at the 
bottom of the band since the same divergence occurs in 
the ladder geometry where (S) recovers to one as p —>■ 0. 

Figures [T](b) and[T](c) show the scaling of the average 
sign with spatial lattice size N and inverse temperature 
P respectively. After plateaux at small N and /3, where 
(S') = 1, the average sign decreases in a manner which is 
largely consistent with exponential. ln(S) is not perfectly 
linear in P or N, but exhibits some downward curvature, 
which we believe indicates the scaling regime has not yet 
been fully attained. We will remark on this more fully 
later in this section. As we shall see, this exponential 
decrease is also the case in other geometries, although the 
decrease with N once one exits the plateau is in general 
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FIG. 1: (Color online) Sign problem for chains of different 
lengths at [/ = 4 (which equals the bandwidth in this 
geometry), (a) The dependence of (S) on the density, (b) 
In (S) vs the system size, N, and (c) vs /3 at several fixed 
densities p. The densities used in panels (b) and (c) are 
indicated by vertical dashed lines in panel (a). The errorbars 
in this figure and in the remaining figures in this paper can 
be inferred from the scatter in the data, and hence, may not 
have been shown. Axis label (sign) is referred to as {S) in 
text. 


less abrupt than with /3. In Fig.[T](b) at p = 0.625,0.875, 
the average sign, despite its exponential decay, remains 
quite manageable out to iV > 200. Even at p = 0.2 the 
spatial size must be tripled from iV ~ 70 before In(S') ^ 
—4. In Fig.jljc), on the other hand, the decay to In(S') ^ 
—4 takes place after only a 50% increase in [3 (from /3 = 8 
to (3 = 12). 

Next, we turn to ladder geometries, which are natural 
extensions of chains, before studying the square lattice. 
Ladders are of interest for several reasons. First, they 
have been extensively studied by DMRG^ as a stepping 
stone to 2D. Second, by changing the ratio of the 
rung hopping to the hopping along the chains, one can 
access U = 0 states that are metallic or band-insulating 
at half-hlling. The effect of these phase changes on the 
sign problem for [/ ^ 0 is one goal of the data presented 
here and in Sec. IVB. 

The results for the average sign on different ladders are 
plotted in Fig. [2] The rung hopping is set to tj. = 1, so 
that the U = 0 band structure is metallic. As with Fig.[l] 
the top panel shows the density dependence for different 





FIG. 2: (Golor online) Same as Fig. [T] but for the ladder 
geometry with the rung hopping equal to the hopping along 
the chains and different values of the interaction strength and 
P, which are noted in the panels. 


system sizes at fixed U,(3 while the bottom two panels 
show that the decay is consistent with exponential in N 
and f3. (See remarks to follow on challenges in capturing 
fully linear behavior of ln(5')(/?).) (S) has a minimum 
at filling p Ri 0.8, similar to what is known to occur also 
for a square lattice. (See Fig. [H) For t±/t > 2 the 
noninteracting system is a band insulator (BI). Because 
the BI occurs at the particle-hole symmetric density p = 
1, (S') is pinned at unity. Therefore, for this case, we will 
examine the signs of the determinants of the individual 
spin matrices (S^) = (S^) in Sec. IVB. 

The square lattice is the most well-studied Hubbard 
model geometry, owing to its possible relevance as a 
simple model of cuprate superconductivity and d-wave 
pairing driven by antiferromagnetic fluctuationsi^ The 
total sign for the square lattice at D = 6,/3 = 4 is 
shown as a function of hlling p for different lattice sizes 
in Fig. [31(a). As mentioned earlier, the peak in the 
4x4 lattice occurs as hve of the 16 allowed k points 
(corresponding to a density p = 10/16 = 0.625) fill up 
prior to half-filling. This peak is even more evidentiS 
at 17 = 4 providing further evidence that, in this case, 
{S) is connected to the shell structure at 17 = 0, though 
the connection appears to diminish with larger 17. It is 
possible that the shell structure would appear on larger 
lattices if lower temperatures were accessible. A rather 
universal feature is the minimum in (S') at p ~ 0.85, 
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FIG. 3: (Color online) Same as Fig. [T] but for the square 
geometry and with different parameters, which are noted in 
the panels. 


which is shared by the ladder geometry. 

The very nature of the sign problem makes it 
challenging to provide a completely compelling linear 
plot of In (S) as a function of p. For the data in 
Fig. 3c, for example, at /3 = 5.33 a run with 10® sweeps 
through the space-time lattice of Hubbard-Stratonovich 
variables (100 spatial sites and 64 time slices) takes 
several hours on a work-station, and gives (S) = 0.0095± 
0.0012 (corresponding to In {S) = —4.66. The slope 
dln(5')/d/3 ^ —4, so we can roughly estimate In (S') ~ 
-8.9 at /3 = 6.33, and hence (S) ~ 0.00013. A 
measurement of this value to 10% accuracy would require 
an error bar of about 0.00001, a factor of ~ 10^ less than 
the error obtained at P = 5.33. Since error bars only go 
down as the square root of the number of measurements, 
such a run would entail ~ 10"*’ times as many sweeps, and 
a cpu time of several months. 

The U dependence of the average sign at fixed /3 = 4 
and N = 100 for a square lattice is shown in Fig. 0] at 
p = 0.625 and 0.875. The evolution of (S) with U is 
rather similar to that with /3: A plateau at weak coupling 
where (S) Ri 1 is followed by an abrupt downturn. Thus, 
in practice, once the sign begins deviating from unity 
there is only a narrow window of stronger couplings where 
data can be acquired. This is demonstrated in Fig. Uas 
(S') as a function of U is shown to decrease exponentially 
once U >A. 



FIG. 4: (Color online) The U dependence of the sign problem 
for the square lattice. 


Here, we briefly generalize our ladder results by 
providing a more complete description of the behavior 
of (S) as the aspect ratio of the lattice goes from one to 
two dimensions. Fig.[S]shows data at 17 = 6 and /3 = 6 for 
lattices with a fixed N = 96 but different aspect ratios. 
As mentioned earlier, the chain geometry seems to be 
rather unique. For all other cases, the average sign is 
close to unity for a range of low densities, and then steps 
downwards to small values in a region centered at p ~ 0.8, 
recovering only at the PHS point p = 1. 

The precipitous nature of the decrease in (S) near 
p = 1, evident in Fig. [5l is quantified in Fig. [6l Panel 
[HKa) shows the logarithm of the sign versus the doping 
away from half-filling. The linear behavior indicates that 
(S) ~ The decay constant a is large and negative. 

Its P and U dependences are given in Fig. |6l(b). To our 
knowledge, this behavior of {S) had not been studied 
before. However, scaling forms for physical observables 
like the compressibility k = dp/dp, as one exits the 
Mott phase at p = 1 have been suggested .®**’®^ In these 
theories, k follows a power law koc(1 — p)“'',so that it 
diverges just before it vanishes. This also occurs in the 
boson Hubbard model.— 

We conclude this section on hypercubic geometries by 
showing the behavior of (S') for a cubic lattice in Fig. [3 
Because the number of lattice sites grows so rapidly with 
linear size, we consider cases where ^ Ly ^ (while 
keeping all linear lengths even to avoid frustration of 
antiferromagnetic correlations.) The qualitative behavior 
is almost identical to that of rectangular lattices, with 
a deep minimum in (S) upon doping from half filling, 
followed by a recovery at p < 0.6. Curves for different 
sizes almost coincide for these large dopings, indicating a 
very slow decay with TV, as seen in Fig. [7Kb). Indeed, 
for some densities {S) even increases as TV increases. 
Presumably, this is a transient phenomenon associated 
with the rather small linear lattice lengths which are 
accessible in 3D. The decay with /3, Fig.jTKc), is, as usual, 
































FIG. 5: (Color online) The average sign as a function of p for 
the rectangular lattices with the same lattice size N = 96 but 
different lengths of the sides. 

rapid. 


B. Other Bipartite Lattices 

Hypercubic lattices are just one instance of bipartite 
geometries which are free of the sign problem at half¬ 
filling owing to the PHS. Here we present data for 
three additional bipartite geometries, all of which are of 
interest because of their materials applications. 

We consider first the “Lieb lattice” which consists of an 
underlying square array of A sites with additional two¬ 
fold coordinated B sites on each bond (see the left panel 
of Fig. [5]). This structure is a more chemically realistic 
multiband representation of the Cu 02 planes of the 
cuprate superconductors, with the Cu atoms forming the 
square array and bridging O atoms. The relevant filling 
of such a three-band model for the cuprates consists of 
one hole per Cu 02 unit cell, that is, well away from half¬ 
filling. Indeed, a realistic model of the cuprates would 



FIG. 6: (Color online) (a) By doping the system away from 
half filling on the square lattice, the average sign decreases 
exponentially, i.e., (S) ~ In (b) the decay constant, 

a, is plotted as a function of U for different values of /3. The 
lattice size is = 10 x 10. 





FIG. 7: (Color online) Same as Fig. [T] but for the cubic 
geometry and with different parameters, which are noted in 
the panels. 



FIG. 8: (Color online) The Lieb lattice (left) and the 1/5- 
depfeted square lattice (right). The latter is taken from 
Ref. [^. Arrows (dashed squares) show the unit vectors (cells) 
for each geometry. 


incorporate a substantial energy epu which represents 
the additional cost of holes to occupy an O p-orbital 
compared to a Cu d-orbital. 

Nevertheless, the half-filled case, with three fermions 
per unit cell and Cpd = 0, has considerable interest: Lieb 
has givenS a rigorous demonstration of a ferrimagnetic 
ground state in this situation. The crucial observation 
is that the numbers of sites on the A and B sublattices 
(iV _4 and TVb, respectively) are unequal. Lieb showed 
that for any bipartite lattice with > N_a there is a 
‘flat band’ with A/g — zero energy eigenstates^. A 
recent DQMC study has explored the attractive Hubbard 
model in this geometryS^. The presence of a flat 











































band was found to have important effects on physical 
properties like the local moment and pairing correlations. 
Figure|9]examines the sign problem for the repulsive case. 
No qualitative difference is discernible from the square 
lattice. Indeed, the sign shows no signal whatsoever as 
it passes through p = 2/3, the filling which corresponds 
to entry into the fiat band. 

The honeycomb lattice is another bipartite lattice we 
study here. It has an interesting semi-metallic density 
of states which vanishes linearly at w —>■ 0. Like the 
hypercubic lattices, it has Nj\^ = Ns- Figure [TUI exhibits 
the average sign in the usual array of panels. (S) is a 
bit reduced in densities p ^ 0.6 compared to the other 
bipartite lattices, but otherwise behaves in a manner 
rather similar to them. 

Our final bipartite geometry is a 1/5-depleted square 
lattice. This is a cousin of the Lieb lattice, in that it can 
be regarded as a square lattice with 1/4 (rather than 1/3) 
of the sites removed, and is the geometry appropriate to 
the magnetic V atoms in CaV 409 (see the right panel of 
Fig. [8]). As with the Hubbard model on a Lieb lattice, 
this model exhibits interesting magnetic orderings. In 
particular, at half-filling and U = 6, as the ratio t'/t of 
the inter- to intra-plaquette hopping is increased, one 
goes from a plaquette singlet phase to a phase with 
antiferromagnetic long-range order at (t'/Oci ~ 0.7 and 
then to a dimer singlet phase where long-range order 
is again absent at {f /t)c 2 ~ 1-3. Figure [TT] shows 





FIG. 10: (Color online) Same as Fig.[T]but for the honeycomb 
geometry and with different parameters, which are noted in 
the panels. 







FIG. 9: (Color online) Same as Fig. [T]but for the Lieb lattice 
and with different parameters, which are noted in the panels. 


the doping dependence of the average sign for the 1/5- 
depleted square lattice for t' = t, which corresponds to 
the ordered phase at p = 1. 

Figure [12] gives the dependence on t'/t for two fixed 
chemical potentials which correspond to p Ri 0.9 and 
0-6 ^ P ^ 0.7. Although the density is varying a bit with 
t'/t, there is a steady decrease in {S) as t'/t decreases. 
The dimer singlet phase at large t'/t has a well-behaved 
sign, while the plaquette singlet phase has a much smaller 
average sign, presumably as a consequence of the fact 
that the sign problem of an isolated 2x2 plaquette is 
much worse than a dimer. 


C. Triangular and Kagome Lattices 

We conclude our survey of lattice geometries with two 
non-bipartite cases: the triangular and Kagome lattices. 
Since there is no PHS, we expect there will be a sign 
problem at all densities, including half filling. Moreover, 
we must present data over a full range of fillings, 0 < 
p < 2. The sign of the hopping is also relevant for these 
structures. Our choices are t = 1 for the triangular lattice 
and t = —1 for the Kagome lattice. In terms of the 
density of states, these choices mean that the DOS is 
nonzero in the ranges [-6, 3] and [-2, 4] for the triangular 
and Kagome geometries, respectively. 

Figure [131 shows results for the triangular lattice. Shell 
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FIG. 11: (Color online) Same as Fig. [T] but for the 1/5- 
depleted square lattice with equal inter- and intra-plaquette 
hopping amplitudes and different parameters for U, j3 and N, 
which are noted in the panels. 


FIG. 13: (Color online) Same as Fig. [T]but for the triangular 
geometry and with different parameters, which are noted in 
the panels. 


structure is evident for the smallest cluster (iV = 9), 
however, as with the square lattice, is absent for larger 
lattices {N = 36,81). The most marked difference from 
the bipartite cases is that (S) is not pinned at one for p = 
1, but otherwise the behavior of {S) is quite similar to the 



FIG. 12: (Color online) Evolution of the average sign on the 
1/5-depleted square lattice as the ratio of the inter- and intra- 
plaquette hopping amplitudes {t'/t) increases. The results are 
obtained for two different fixed chemical potentials, p — —1.5 
and -2.5. The filled symbols show the average sign and the 
empty symbols denote the evolution of the average density 
p. Note that p varies slowly in most of the range of t'/t and 
only increases significantly towards one for t/t' < 0.5. The 
bandwidth is kept fixed at 6 for all values of t'/t. 


previous cases. The same is true of the Kagome lattice 
in Fig. 1141 except that there is a persistent bump in {S) 
for p slightly less than quarter filling. (A similar feature 
was noted for chains.) The structure in the average sign 
is somewhat more asymmetric about half filling than the 
triangular case. One feature which does not seem to be 
shared with other geometries is the existence of an abrupt 
change in (S') appearing here at p ~ 1.1, so that it is 
well behaved for most densities p > 1.1. The scalings 
with the cluster size and [3 of the average sign for these 
non-bipartite geometries, shown in the bottom panels of 
Figs. [13] and [Ml follow a similar pattern as for the other 
bipartite geometries. 


IV. FURTHER ANALYSIS 


We now discuss possible patterns which emerge from 
these datasets. We focus on three areas: the role 
of the density of states, the contribution of individual 
spin components to the sign problem, and spatial 
entanglement. A fourth feature of (S'), scaling in the 
vicinity of the PHS point, was discussed previously. 
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FIG. 14: (Color online) Same as Fig. [T] but for the Kagome 
geometry and with different parameters, which are noted in 
the panels. 

A. Role of the U = 0 Density of States 

It seems plausible that the noninteracting density of 
states could play an important role in the sign problem. 
In this subsection we make a few observations on that 
possibility. 

The square and honeycomb lattices have quite 
dramatically different U = 0 densities of states, especially 
near w = 0 where the DOS diverges logarithmically 
for the square lattices, and vanishes linearly for the 
honeycomb lattice. Yet, if we compare the behaviors of 
(S) as a function of filling in Figs. 15] and [TUI we see little 
qualitative difference. Both evolutions exhibit a rapid 
fall-off from (S') = 1 at the PHS point p = 1, a broad 
minimum centered at p ^ 0.8, followed by a recovery to 
(S) = 1 in the dilute limit. The differences in DOS are 
even more diverse among the other geometries studied 
here. However, the special features of the DOS, which 
could trend in completely opposite directions, appear to 
have little to no effect on the behavior of (S) here. 

Indeed, we have noted already that, more generally, 
the different geometries and their associated distinct 
densities of states all share a qualitatively similar 
behavior of (S) with doping. The only ‘unique’ geometry 
was the one-dimensional case where, for example (S) did 
not recover to one at low densities. 

Most of the data we present are for interaction 
strengths U at least four times the fermion hopping t. 



FIG. 15: The product of the expectation values of the signs 
of the individual spin np and spin down matrices is shown 
at half-filling as a function of the ratio of rung to on-chain 
hopping for the ladder geometry. ( (S-|-) = {Si) by symmetry.) 
The signs are well-behaved in the band insulating phase 
t±/t > 2. 


and hence, roughly speaking, at least half the bandwidth 
W. These sorts of interaction strengths are the ones 
typically studied in examining magnetic, pair, and charge 
correlations in Hubbard models. The conclusion of the 
observations above seems to be that U > W/2 brings 
the system far enough from the U = 0 limit that 
most features of the noninteracting DOS are no longer a 
controlling factors in the sign problem. Apparently, the 
space and imaginary time fluctuations of the Hubbard- 
Stratonovich held whose effect on the fermions 

increases with U through the parameter A of Eq. 21 
smear out effects of the U = 0 energy levels on (S). 
This appears consistent with a comparison of Fig. EKa) of 
this paper with Fig. 10(a) in Ref. [13 for 4 x 4 Hubbard 
lattices. With (3U constant, the sharpness of the feature 
near p = 0.6 is signihcantly reduced when U is increased 
from 4 to 6. 

This lack of dependence on the DOS is the case even 
for a hat band, where the very large delta function in 
the DOS might have been expected to have an especially 
discernible impact on {S). However, in the case of the 
Lieb lattice (Fig. 0, (S) behaves completely smoothly 
through the edge of the hat band at p = 2/3. We do 
note that broad features in the U = 0 DOS do appear to 
have some correlation with the behavior of (S^) = (Si). 
For example, for the chain, square, and cubic lattices, we 
observe that (Sa) tends to be smaller when a “smoothed” 
DOS is larger (not shown). Behavior near half-hlling for 
bipartite lattices is additionally mediated, as discussed 
previously, by the fact that at half-hlling (S') = 1 by 
symmetry. 
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FIG. 16: (Color online) The product of the average signs for 
each spin specie on the 1/5-depleted square lattice at half 
filling as a function of t' jt. 


B. Spin Components 

Thus far, we have focused almost exclusively on the 
total sign (5). As commented in Sec. IIB, assuming there 
is no ‘off-diagonal’ term in the Hamiltonian, allowing 
the mixing of different spin species, the fermionic trace 
results in separate determinants for each a. Here, we 
do a further analysis of the signs of the individual spin 
components. We have two goals: first, we would like 
to examine the correlations between the up and down 
determinants for the same HS configurations. Second, 
as in the case of the ladder or the 1/5-depleted square 
geometry, the spin-resolved sign may offer some insight 
into the potential connection of the average sign to other 
observables and phase transitions that occur in the PHS 
regimes (half filling). 

We start with the second objective. We consider a 
ladder geometry with intra-chain hopping t and inter¬ 
chain (rung) hopping t±. The U = 0 band structure is 
e±{kx) = ±t_L — 2tcoskx, so that the ladder is a HI for 
t±/t > 2 at half-filling, and a metal for t±/t < 2. Figure 
[m shows the product {S^){Si) as a function of t±/t for 
p = 1 (where (S^Si) = 1) and U = 4 for /3 = 4, 5 and 
6 . There seems to be some evidence that entering into 
the metallic phase at t±/t = 2 coincides with an increase 
in the number of negative determinants. The t±/t = 0 
limit is also interesting. It corresponds to two deconpled 
ID chains. Evidently, even the individual determinants 
are free of negative signs at half-filling. 

Unlike for the ladder geometry, the product of the 
average signs for the two spin species does not display 
any signature at the (magentic) phase transitions for the 
1/5-depleted square lattice. This can be inferred from 
Fig.lTClwhere {S^){Si) at half filling is plotted vs t'/t for 
U = 6 and /3 = 4. Similar to the total sign away from 
half filling (Fig. [12]), this quantity also decreases from one 
in the large t'/t decoupled dimers region to near zero in 
the small t' jt plaquette region. There are two quantum 
phase transitions^^ for this value of the interaction at 
t'/t ~ 0.65 and t/t' ^ 0.77. However, the average signs 
vary smoothly around these values with no special feature 
(at this temperature) that can be attributed to the phase 




FIG. 17: Correlations between the average of the product 
of the up and down determinant signs and the product of 
the averages. Top panel: square lattice (PHS case). Bottom 
panel: triangular lattice (non-PHS case). The data indicate 
that, for the most part, the signs of the up and down 
determinants are independent, except at the PHS density 

p = 1. 


transitions in the ground state. 

With this in mind, we analyze the correlations between 
the up and down average signs for the square and 
triangular lattices in Fig. |T7| In the former geometry, at 
half-filling, because of PHS there is a perfect correlation 
of the up and down signs, i.e., (S') = 1, which differs 
substantially from Correlations remain in 

the neighborhood of the PHS point. However, for 
densities p < 0.8 (and, by symmetry, p > 1.2), (S) = 
(S'|-)(S^) to a high degree of accuracy, indicating that 
the signs of the spin up and spin down matrices for 
given HS configurations are essentially independent. It 
is interesting that the density at which correlations 
disappear roughly coincides to where the sign problem 
is the worst. 
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We expect the absence of the PHS to alter the 
relationship between the total sign and its individual 
components at half filling. This issue is examined in the 
triangular lattice results in the bottom panel of Fig. [T71 
The determinant signs are generally independent over 
the entire density range, with {S) = {S^){Si), except 
for a small region around three quarter filling, p fa 1.5. 
This is evidently related to the bumps in (S) already 
noted in Fig. 1131 and provides additional insight into 
that phenomenon: The increase in (S') is associated, 
at least in part, with a larger correlation of the two 
determinant signs. Interestingly, the latter seems to be 
taking place close to where the peak in the DOS is, 
suggesting a possible connection, however, a similar but 
less prominent bump also appears around p ^ 0.7, where 
there is no feature in the DOS. 


C. Spatial Entanglement 

Our data have supported the exponential decay of the 
sign with N (Fig. [3 for the cubic lattice is inconclusive, 
perhaps due to the shorter lengths of the sides). Here, we 
point out that there is a limit in which this phenomena 
can be rigorously established. Specifically, we observe 
that the sign must decay exponentially with N for lattices 
consisting of decoupled clusters (by, for example, setting 
the inter-cluster hopping to zero). Clearly the average 
sign of the entire system will be the product of the signs 
of the subclusters. If the subclusters are identical, an 
exponential decrease of the average sign with system 
size (number of subclusters) trivially follows. One can 
then ask what is the effect of linking the clusters. If 
the connection is through non-zero inter-cluster hopping, 
a ‘world-line’ picture would generally suggest the sign 
should get even smaller, as the fermion paths are allowed 
additional opportunity to exchange. In fact, precisely 
this would happen as a lattice of decoupled chains in the 
X direction with zero hopping in the y direction, ty = 0, 
is converted into a 2D lattice by turning on ty. The 
ty = 0 limit has no significant sign problem in world-line 
approaches, whereas large ty does. 

In DQMC, the opposite can occur: the linking 
(‘entanglement’) of spatial clusters mitigates the sign 
problem. This can be seen in several Hubbard 
Hamiltonian geometries which contain decoupled clusters 
as a limiting case. The simplest is the ladder 
geometry where t_L = 0 corresponds to two independent 
chains. Other examples are the ‘plaquette’ model 
studied by Scalapino^, Kivelson^ and otherS)^^— 
as a description of superconductivity arising from pair 
binding on 2 x 2 plaquettes, and the 1/5-depleted square 
lattic o^^i^^i^^ for which the results in the left panel of 
Fig. [T3 show the improvement of the sign by increasing 
the inter-plaquette hopping t' from zero. 


V. CONCLUSIONS 

The sign problem remains one of the fundamental 
challenges in computational physics. This paper 
has focused on the sign problem in determinant 
quantum Monte Carlo. Our goal has been to bring 
together data for a collection of geometries (hypercubic, 
ladder, depleted square, Lieb, honeycomb, Kagome, 
and triangular lattices) and parameter (temperature, 
interaction strength, and density) ranges that are, at 
present either not available or, at best, scattered through 
the literature. There are other cases one could study. 
However, this extensive set already enables us to make 
some general inferences about the sign problem, at least 
in the specific case of the DQMC, which may have other 
applications as well. 

We first considered the general behavior of the sign (S') 
as a function of inverse temperature /3 = l/^ksT) and 
of lattice size N. Arguments for exponential behavior 
have been made in this regard involving the winding of 
world lines, in approaches where such paths are sampled 
in the simulation. However, it has been noted that this 
reasoning does not necessarily transfer to the auxiliary 
field approach of DQMC, where it is difficult to see 
how to map the problem onto a sum over world line 
paths of an effective Hamiltonian with local interactions. 
Without such a theoretical basis, it is of interest to 
explore numerically the scaling of (S') with /? and N for 
a variety of cases, to determine the general behavior. 

Previous work regarding /3 used the Hubbard model 
on 3 X 2 lattices^ and 4x4 lattices^ii finding an 
exponential decay of (S) with /3 for large /3. Our data, 
on considerably larger spatial lattices and for a variety 
of different geometries, is also consistent with a large 
(5 exponential decay, with a small /3 regime throughout 
which (S) ~ 1, and eliminates ‘shell effects’ seen in 
smaller lattices. 

We next explored the scaling of (S) with lattice size 
N. Previous work, comparing 2D 4 x 4, 6 x 6, and 8x8 
Hubbard lattices^d^ and 3D 4x4x4 and 6x6x6 
lattices^ was inconclusive. In both the 2D and 3D 
cases, for example, for some parameter ranges (S') was 
actually worst for the smallest systems studied. Using 
larger lattice sizes and differing geometries, our results 
are clearly consistent with an exponential decay of (S) 
with N for large N. In a few but not all cases, we also 
found a small N regime in which (S) fa 1, similar to what 
was seen for j3. Data based on an examination of (S) as a 
function of U were consistent with a similar exponential 
decay of (S') with U, for U sufficiently large. 

As is well known, the sign problem is ‘cured’ in special 
particle-hole symmetric cases like the half-filled Hubbard 
model on a bipartite lattice, where the determinants 
always come in pairs which share the same sign. This 
paper presented a more general study of correlations 
between the signs. Near half-filling, for the square lattice, 
we found (S) > (S-j-) (Sj,), where (Sa) is the average sign 
of a single spin determinant, see Eq. 6. This indicates 
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the development of strong correlations between the signs 
of the up and down determinants, 

A closely related issue is the precise form of the decay 
of (S) close to the half-filled p = 1 PHS point where 
(S) = 1. We found a rapid decay, consistent with the 
form (S) = and that “a”becomes increasingly 

large and negative with U and with /3. This indicates 
a rapid loss of the positive correlations between up and 
down determinants as one moves away from half filling. 

Lastly, we considered the effect of “entanglement”, 
where initially decoupled clusters were linked together 
with fermion hopping terms. This linking increased {S) 
improving the sign problem, in complete contrast to what 
one would expect in a “world line” picture, where the 
additional windings available due to the extra hopping 
terms would typically be expected to worsen the sign 
problem. 

Another significant feature of this work is that 
the geometries being studied have a wide range of 
noninteracting densities of states. A crucial conclusion of 
our work, evident in comparing the top panels of Figs. 2- 
7 is that, despite the wide variation in the DOS, on 
bipartite lattices, the behavior of the sign with density 
is almost “universal” in the sense that it always falls 
rapidly away as the lattice is doped, attains a minimum 
in the vicinity of p ~ 0.8 and then recovers. This 
observation is of interest because of the rather different 
physics expected of strongly correlated electrons on these 
geometries. The only exception to this universality is 
seen in Fig. [U for the one dimensional chain. 

Overall, features in the U = 0 DOS associated with 
special behavior at an isolated energy seem to have little 
effect on (S). This may be due to the “smoothing out” of 
these features by the HS field fluctuations which increase 
with increasing U. There does, however, seem to be some 
indication that a larger value of a “smoothed” DOS is 
associated with smaller average signs (S^) = (S'j,) of the 
determinants of the matrices of individual spins. 

Another interesting possibility for further exploration 
is that of a linkage between (S) and the spectral 
function at the Fermi level, A{uj = 0). Here A(uj) 
is related to the time dependent Green’s function via 
G(r) = /dr A{lo) / (-b 1) and equals the 
density of states N(lo) in the non-interacting limit. The 
spectral function incorporates the effects of U and hence, 
potentially, might correlate better with the sign. There 
are some hints that this is the case, for example, on a 
square lattice the U = 0 density of states diverges at half¬ 
filling yet the low-temperature Aioj) vanishes for all U 


and the sign behaves perfectly. Indeed, the sign behaves 
well in a range of chemical potentials p within the Mott- 
Slater gap. One could imagine pursuing this possible 
connection more closely through computing the spectral 
functions for the various geometries and parameter values 
considered here. However, this would be a major 
task beyond the scope of this paper, involving the 
measurement of the imaginary-time dependent Green’s 
function and its analytic continuation to real space. In 
addition, because it is more difficult to obtain A^uj), as 
compared to the noninteracting density of states, it is 
unclear what the utility of the discovery of a connection 
between A{uj) and the sign would be. 

Despite the data and interpretation presented here, 
the sign problem remains a big mystery. Of particular 
interest is the possible relation between {S) and the 
underlying physics of correlated electrons. For example, 
early in the development of DQMC, Hirsch pointed out^i 
a mapping between the spin-spin correlation function, a 
property solely of the Hubbard Hamiltonian itself, and a 
correlation function of the Hubbard-Stratonovich field X. 
Since the fluctuations in X ultimately determine the sign 
problem, this suggests the possibility that the behavior 
of {S) might be related to some appropriate observable. 

We end on a speculative note in this regard. The 
unfortunate coincidence (if it is a coincidence) that for 
the square lattice (S') is worst behaved very close to the 
most interesting ‘optimal’ doping, where is largest 
in the cuprate superconductors, has often been opined. 
One thing we observe here is that this dip in (S) at 
p ft! 0.8 is absent in d = 1 but appears already in 
the ladder geometry. Since Hubbard ladders appear 
to show signatures of d-wave pairing^ the possibility 
of a deep connection between the fermion sign and 
superconductivity in the Hubbard Hamiltonian remains 
a possibility. 
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